﻿using MathNet.Numerics.Statistics;
using phminarm.algorithm.Base;
using phminarm.algorithm.DataProcess;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Reflection;
using System.Text;
using System.Threading.Tasks;

namespace phminarm.algorithm.FeatureExtraction
{
    public class FrequencyDomainFeatures
    {
        private float fs;//采样频率
        private float df;//频率分辨率
        private int N;//采样点数
        private float[] timewave;
        private float[] spectrum;
        private float[] phase;
        public FrequencyDomainFeatures(float[] timewave, float fs, WindowType windowType = WindowType.Hanning)
        {
            this.fs = fs;
            this.N = timewave.Length;
            this.df = fs / N;
            this.timewave = timewave;
            (this.spectrum, this.phase) = FFTWraper.GetFFTAmp(timewave, windowType);
        }

        public FrequencyDomainFeatures(float[] spectrum)
        {
            this.timewave = new float[0];
            this.spectrum = spectrum;
        }
        //1.幅值均值
        public float meanamp { get; private set; }
        //2.幅值方差
        public float varinceamp { get; private set; }
        //3.幅值标准差
        public float standardamp { get; private set; }
        //4.幅值偏度因子
        public float skewnessfamp { get; private set; }
        //5.幅值峭度因子
        public float kurtosisfamp { get; private set; }
        //6.重心频率
        public float barycenterfreq { get; private set; }
        //7.均方频率
        public float msfreq { get; private set; }
        //8.均方根频率
        public float rmfreq { get; private set; }
        //9.频率方差
        public float varincefreq { get; private set; }
        //10.频率标准差
        public float standardfreq { get; private set; }
        //11.频率分布系数
        public float distributionfreq { get; private set; }
        //12.频率偏度
        public float skewnessfreq { get; private set; }
        //13.频率峭度
        public float kurtosisfreq { get; private set; }
        //14.平方根比率
        public float mrratio { get; private set; }
        //15.轴承频率
        public List<Bear> bears { get; private set; }
        public void GetBearFeature(List<Bear> bears, float[] rpm)
        {
            this.bears = bears;
            float span = 1;//频率搜索范围±1Hz
            for (int i = 0; i < bears.Count; i++)
            {
                //计算轴承特征频率
                bears[i].GetFeatures(rpm[i]);
                //提取轴承特征频率幅值
                bears[i].BPFOA = findMax(bears[i].BPFO, span);
                bears[i].BPFIA = findMax(bears[i].BPFI, span);
                bears[i].FTFA = findMax(bears[i].FTF, span);
                bears[i].BSFA = findMax(bears[i].BSF, span);
            }
        }
        //16.齿轮频率
        public List<Gear> gears { get; private set; }
        public void GetGearFeature(List<Gear> gears, float rpm)
        {

        }
        //17.转速倍频
        public float[] SpeedMultiplication { get; private set; }
        public void GetSpeedMultiplication(float rpm, int x)
        {
            float[] SpeedMultiplication = new float[x + 1];
            float fr = rpm / 60;
            float fspan = 1;//搜索频率范围为±1Hz
            SpeedMultiplication[0] = findMax(fr * 0.5f, fspan);
            for (int i = 1; i <= x; i++)
            {
                SpeedMultiplication[i] = findMax(fr * i, fspan);
            }
        }
        //18.振动加速度特定频带范围内的RMS
        public float rms_acc_10_10k { get; private set; }
        //19.振动速度特定频带范围内的RMS
        public float rms_vel_10_1k { get; private set; }

        //提取全部频域特征
        public void GetAllFrequencyFeatures(float lowerFreq = 10, float upperFreq = 10000)
        {
            int lower = Convert.ToInt32(Math.Floor(lowerFreq / df));
            int upper = Convert.ToInt32(Math.Ceiling(upperFreq / df));
            int M = upper - lower + 1;
            float[] data = spectrum.Skip(lower).Take(M).ToArray();
            //1.幅值均值
            meanamp = data.Average();
            //2.幅值方差
            varinceamp = Convert.ToSingle(data.Variance());
            //3.幅值标准差
            standardamp = Convert.ToSingle(data.StandardDeviation());
            //4.幅值偏度因子
            skewnessfamp = Convert.ToSingle(data.Sum(x => Math.Pow((x - meanamp), 3)) / (M - 1) / Math.Pow(standardamp, 3));
            //5.幅值峭度因子
            kurtosisfamp = Convert.ToSingle(data.Sum(x => Math.Pow((x - meanamp), 4)) / (M - 1) / Math.Pow(standardamp, 4));
            //6.重心频率
            barycenterfreq = 0;
            for (int i = 0; i < M; i++)
            {
                barycenterfreq += data[i] * (i + lower) * df;
            }
            barycenterfreq = barycenterfreq / data.Sum();
            //7.均方频率
            msfreq = 0;
            for (int i = 0; i < M; i++)
            {
                msfreq += Convert.ToSingle(data[i] * Math.Pow((i + lower) * df, 2));
            }
            msfreq = msfreq / data.Sum();
            //8.均方根频率
            rmfreq = Convert.ToSingle(Math.Sqrt(msfreq));
            //9.频率方差
            varincefreq = 0;
            for (int i = 0; i < M; i++)
            {
                varincefreq += Convert.ToSingle(data[i] * Math.Pow((i + lower) * df - barycenterfreq, 2));
            }
            varincefreq = varincefreq / data.Sum();
            //10.频率标准差
            standardfreq = Convert.ToSingle(Math.Sqrt(varincefreq));
            //11.频率分布系数
            distributionfreq = varincefreq * data.Sum() / (M - 1);
            //12.频率偏度
            skewnessfreq = 0;
            for (int i = 0; i < M; i++)
            {
                varincefreq += Convert.ToSingle(data[i] * Math.Pow((i + lower) * df - barycenterfreq, 3));
            }
            varincefreq = Convert.ToSingle(varincefreq / M / Math.Pow(distributionfreq, 3));
            //13.频率峭度
            kurtosisfreq = 0;
            for (int i = 0; i < M; i++)
            {
                kurtosisfreq += Convert.ToSingle(data[i] * Math.Pow((i + lower) * df - barycenterfreq, 4));
            }
            kurtosisfreq = Convert.ToSingle(kurtosisfreq / M / Math.Pow(distributionfreq, 4));
            //14.平方根比率
            mrratio = 0;
            for (int i = 0; i < M; i++)
            {
                mrratio += Convert.ToSingle(data[i] * Math.Sqrt((i + lower) * df - barycenterfreq));
            }
            mrratio = Convert.ToSingle(mrratio / M / Math.Sqrt(distributionfreq));

            //18.振动加速度特定频带范围内的RMS
            rms_acc_10_10k = 0;
            rms_acc_10_10k = calculaterms(spectrum, df, 10, 10000);
            //19.振动速度特定频带范围内的RMS
            rms_vel_10_1k = 0;
            //频域积分
            float[] vel = Integration.FrequencyDomain(spectrum, df, 0);
            rms_vel_10_1k = calculaterms(vel, df, 10, 1000);
        }

        /// <summary>
        /// 搜索范围内的最大值
        /// </summary>
        /// <param name="x">搜索频率点Hz</param>
        /// <param name="dx">搜索范围Hz</param>
        /// <returns></returns>
        private float findMax(float x, float dx)
        {
            if (x < 0 || dx < 0) return -1f;

            int n = Convert.ToInt32(Math.Round(x / df));
            int dn = Convert.ToInt32(Math.Round(dx / df));

            if (n > this.N) return -1f;
            if (dn == 0) return spectrum[n];

            //skip和take函数已经考虑了数组越界的情况
            float[] span = spectrum.Skip(n - dn).Take(2 * dn).ToArray();
            return span.Max();
        }

        private float calculaterms(float[] spectrum, float df, float lowerFreq, float upperFreq)
        {
            double rms = 0;
            int lowerN = (int)Math.Floor(lowerFreq / df);
            int upperN = (int)Math.Ceiling(upperFreq / df);

            for (int i = lowerN; i <= upperN; i++)
            {
                rms += Math.Pow(spectrum[i] / FFTWraper.recover * 1.633 / Math.Sqrt(2), 2);
            }

            rms = Math.Sqrt(rms);

            return (float)rms;
        }
    }
}
